/*
*		
*		A compilacao e feita por: gcc -o gauss_gmp gauss_gmp.c -lgmp -lm 
*
*/

#include <stdio.h>
#include <stdlib.h>
#include <limits.h>
#include <time.h>
#include <gmp.h> 

#define N 25
#define BASE 10
#define MANTISSA 1024*1024*32


// Inicializa o valor de b, onde a variavel b e passada por parametro
// b = 1/sqrt(2)
void init_mpf_b(mpf_t rop)
{
	mpf_set_default_prec(MANTISSA);
	mpf_t x, y;

	mpf_init(x);
	mpf_init(y);
	
	// x = 1
	mpf_set_ui(x, 1);
	// y = sqtr(2)
	mpf_sqrt_ui(y, 2);
	// x = 1/sqrt(2)
	mpf_div(x, x, y);

	// return x
	mpf_set(rop, x);

	//libera memoria
	mpf_clear(x);
	mpf_clear(y);
}


// calcula a media simples entre a e b, onde as variaveis a,b e 'rop' sao passadas por parametro
// (a + b) / 2
void media_simples(mpf_t rop, mpf_t a, mpf_t b)
{
	mpf_t aux;	
	mpf_init(aux);

	//aux = a + b;
	mpf_add(aux, a, b);
	
	//rop = aux/2
	mpf_div_ui(rop, aux, 2);

	//libera memoria
	mpf_clear(aux);
}


// calcula a media simples entre a e b, onde as variaveis a,b e 'rop' sao passadas por parametro
// sqrt(a + b)
void media_geometrica(mpf_t rop, mpf_t a, mpf_t b)
{	
	mpf_t aux;
	mpf_init(aux);

	//aux = a + b;
	mpf_mul(aux, a, b);
	
	//rop = sqrt(aux)
	mpf_sqrt(rop, aux);

	//libera memoria
	mpf_clear(aux);
}

// calcula o valor de t onde as variaveis t, x, y e a sao passadas por parametro
// t = t - x * (y-a)^2;
void equacao_t(mpf_t t, mpf_t x, mpf_t y, mpf_t a)
{
	mpf_t aux;
	mpf_init(aux);
	

	// aux = y - a
	mpf_sub(aux, y, a);
	// aux = (aux)^2
	mpf_pow_ui(aux, aux, 2);
	// aux = aux * x
	mpf_mul(aux, aux, x);
	//aux = t - aux
	mpf_sub(t, t, aux);

	//libera memoria
	mpf_clear(aux);
}


// calcula o valor de pi onde as variaveis pi, a, b e t sao passadas por parametro
// pi = ((a + b)^2)/(4.0*t)
void equacao_pi(mpf_t pi, mpf_t a, mpf_t b, mpf_t t)
{
	mpf_t aux;
	mpf_init(aux);
		
	// aux = a + b
	mpf_add(aux, a, b);
	// aux = aux^2
	mpf_pow_ui(aux, aux, 2);
	// aux = aux / t
	mpf_div(aux, aux, t);
	// pi = aux/4
	mpf_div_ui(pi, aux, 4);

	//libera memoria
	mpf_clear(aux);
}


int main(int argc, char* argv[])
{
	unsigned long int i; // variavel auxiliar
	time_t t_zero, t_final;
	double delta_t;
	
	// Marca o tempo
	time(&t_zero);
	
	// Float com precisao arbitraria
	mpf_set_default_prec(MANTISSA);
	mpf_t a, b, t, x, y, pi;

	//inicializa as variaveis
	mpf_init(a);
	mpf_init(b);
	mpf_init(t);
	mpf_init(x);
	mpf_init(y);
	mpf_init(pi);

	
	//	a = 1.0;
	//	b = (1.0/(sqrt(2)));
	//	t = (1.0/4.0);
	//	x = 1.0;
	//	y = 1.0; // Independente
	mpf_set_d(t, (1.0/4.0));
	mpf_set_ui(a, 1);
	mpf_set_ui(x, 1);
	init_mpf_b(b);

	
	/*
	 * 	A(n+1) = A(n)+B(n)/2
	 * 
	 *  B(n+1) = sqrt(A(n)+B(n))
	 * 
	 *	T(n+1) = T(n) - P(n)(A(n) - A(n+1))^2
	 *  
	 * 	P(n+1) = 2*P(n)
	 * 
	 * 	pi = ((A(n)+B(n))^2)/4T(n)
	 * 
	 * */

	for(i = 0; i < N; i++)
	{
			mpf_set(y, a);				//y = a;
			media_simples(a, a, b);  	//a = (a + b)/2.0;
			media_geometrica(b, b, y);	//b = sqrt(b * y);
			equacao_t(t, x, y, a);		//t = t - x * pow((y - a),2);
			mpf_mul_ui(x, x, 2);		//x *= 2;
			equacao_pi(pi, a, b, t);	//pi = pow(( a+ b),2)/(4.0 * t);
	}

	//libera memoria
	mpf_clear(a);
	mpf_clear(b);
	mpf_clear(y);
	mpf_clear(t);

	//mostra o valor de pi
	mpf_out_str(NULL, BASE, 0, pi);

	//libera memoria
	mpf_clear(pi);
	
	// mostra o tempo gasto
	time(&t_final);
	delta_t = difftime(t_final, t_zero);
	printf("\nTempo = %.2lf segundos\n", delta_t);

	return 0;
}

